
// TableD2: Heterogeneous effects. Regression coefficients. Couples, 2015-2017.

*===============================================================================
// Step 1: Create dataset - reduce to couples and create heterogeneity variables
*===============================================================================
use    "${dataout}MainDataset" , clear
keep if year>=2015 & year<=2017
keep if couple==1

// The file "Link_grkrets_regtype1" maps basic statistical units to aggregated region types.
// Region types are merged to basic statistical unit of residence.
sort grkrets
merge m:1 grkrets using  "${datain}Link_grkrets_regtype1", keep(match master)
drop _merge

gen CityCategory=. // Creating four categories based on region types
replace CityCategory=1 if (regtype1==1 | regtype1==3) & !missing(regtype1) // Oslo, Bergen, Trondheim or Stavanger
replace CityCategory=2 if (regtype1==2 | regtype1==4) & !missing(regtype1) // Suburbs of the four largest cities
replace CityCategory=3 if (regtype1==5 | regtype1==6) & !missing(regtype1) // Other, smaller cities
replace CityCategory=4 if (regtype1==7) & !missing(regtype1) // Rural areas

// Merging region types to basic statistical unit of the workplace.
forvalues i=1/2 {
	preserve
		use "${datain}Link_grkrets_regtype1" , clear
		rename grkrets grk_bed`i'
		rename regtype1 regtype1_bed`i'
		tempfile helpfile_`i'
		save `helpfile_`i''
	restore
	sort grk_bed`i'
	merge m:1 grk_bed`i' using  `helpfile_`i'', keep(match master)
	drop _merge
}

forvalues i=1/2 {
	gen CityCategory_bed`i' = . // Creating four categories based on region types (workplace of each spouse)
	replace CityCategory_bed`i'=1 if (regtype1_bed`i'==1 | regtype1_bed`i'==3) & !missing(regtype1_bed`i')
	replace CityCategory_bed`i'=2 if (regtype1_bed`i'==2 | regtype1_bed`i'==4) & !missing(regtype1_bed`i')
	replace CityCategory_bed`i'=3 if (regtype1_bed`i'==5 | regtype1_bed`i'==6) & !missing(regtype1_bed`i')
	replace CityCategory_bed`i'=4 if (regtype1_bed`i'==7) & !missing(regtype1_bed`i')
}

* === Interaction variable for geography =======================================
gen commuting_dir = .
replace commuting_dir = 1 if CityCategory == 1 													 // Within city
replace commuting_dir = 2 if CityCategory == 1 & CityCategory_bed1 != 1 & CityCategory_bed1 != . // From city
replace commuting_dir = 2 if CityCategory == 1 & CityCategory_bed2 != 1 & CityCategory_bed2 != . // From city
replace commuting_dir = 3 if CityCategory == 2 													 // Within suburb
replace commuting_dir = 4 if CityCategory == 2 & CityCategory_bed1 == 1 						 // Suburb to city
replace commuting_dir = 4 if CityCategory == 2 & CityCategory_bed2 == 1 						 // Suburb to city
replace commuting_dir = 5 if CityCategory == 3 													 // Other cities
replace commuting_dir = 6 if CityCategory == 4 													 // Rural areas

* === Interaction variable for income ==========================================
egen wies_fam_sum =rowtotal(wies1 wies2), missing
xtile wies_fam_sum_quintile =  wies_fam_sum, nq(5)
compress

*===============================================================================
// Step 2: Run regressions, four different outcomes, save as ster-files
*===============================================================================
/*
NB: In the do-file "Figure3", the same regression is run for "BEV_fam_yes"
    and "car_fam_count". Therefore, we don't need to run them again.
*/

local outcomelist /*BEV_fam_yes car_fam_count */ car_fam_yes ICE_fam_count 
local treatvar1list toll_fam_mean         
local treatvar2list ptl_fam_km_mean    

** Regressions - same regression new baseline
foreach outcome in `outcomelist'{
	foreach treatvar1 in `treatvar1list'{
		foreach treatvar2 in `treatvar2list'{
			forvalues i = 1/5 { // interactions betw treatment var and income
				gen tr1_inc_q`i' = `treatvar1' * (wies_fam_sum_quintile == `i')
				gen tr2_inc_q`i' = `treatvar2' * (wies_fam_sum_quintile == `i')
			}

			forvalues i = 1/6 { // interactions betw treatment var and geography
				gen tr1_geo`i' = `treatvar1' * (commuting_dir == `i')	
				gen tr2_geo`i' = `treatvar2' * (commuting_dir == `i')	
			}
			compress	
			reghdfe `outcome' 									/// Outcome variable
				`treatvar1' `treatvar2' 						/// Treatvars, baseline values: income quantile 5 and living+working in large city
				tr1_inc_q1 tr1_inc_q2 tr1_inc_q3 tr1_inc_q4  	/// Treatvar1 by income groups
				tr2_inc_q1 tr2_inc_q2 tr2_inc_q3 tr2_inc_q4  	/// Treatvar2 by income groups
				tr1_geo2 tr1_geo3 tr1_geo4 tr1_geo5 tr1_geo6 	/// Treatvar1 by geography
				tr2_geo2 tr2_geo3 tr2_geo4 tr2_geo5 tr2_geo6 	/// Treatvar2 by geography  
				$age  $distance  $time  $publictime  $publicquality, /// Controls
				absorb($FE   $household  $income  $employment  $education i.year) /// Fixed effects and additional absorbed controls
				vce(cluster $clustervar) 						/// multi-dimensional clustering
			summarize `outcome' if e(sample)==1
			estadd scalar MeanDep = r(mean), replace
			summarize `treatvar1' if e(sample)==1
			estadd scalar Treat1 = r(mean), replace
			summarize `treatvar2' if e(sample)==1
			estadd scalar Treat2 = r(mean), replace
			estadd local year="2015-17" , replace
			estadd local grkFE="\checkmark", replace
			estadd local grkbFE="\checkmark", replace
			estadd local HHcon="\checkmark", replace
			estadd local HHWork="\checkmark", replace
			estadd local PublicTime="\checkmark", replace
			estadd local grkFE_W="Yes", replace
			estadd local grkbFE_W="Yes", replace
			estadd local HHcon_W="Yes", replace
			estadd local HHWork_W="Yes", replace
			estadd local PublicTime_W="Yes", replace
			eststo reg_pooled_`outcome'_`y' 
			estimates save "${ster}hetero_reg_O_`outcome'" , replace

			drop tr1_*
			drop tr2_*
		}
	}
}

*===============================================================================
// Step 3: Load ster files, make table
*===============================================================================

eststo clear
estimates use "${ster}hetero_reg_O_BEV_fam_yes"
eststo 
estimates use "${ster}hetero_reg_O_car_fam_count"
eststo 
estimates use "${ster}hetero_reg_O_car_fam_yes"
eststo 
estimates use "${ster}hetero_reg_O_ICE_fam_count"
eststo 

* This table contains the same coefficients as in Table D.2
* However, table D.2 in the online appendix is pivoted manually such that
* - Baseline effect, income groups and geographic areas are columns
* - The four outcome variables are four row-wise panels
* - Each panel contains one row for Road toll (NOK) and one row for bus lane (km)

esttab *, keep(toll_fam_mean /// Baseline toll effect
	tr1_inc_q1 tr1_inc_q2 tr1_inc_q3 tr1_inc_q4 /// Toll, income groups 1-4
	tr1_geo2 tr1_geo3 tr1_geo4 tr1_geo5 tr1_geo6 /// Toll, geography 2-6
	ptl_fam_km_mean /// Baseline bus lane effect
	tr2_inc_q1 tr2_inc_q2 tr2_inc_q3 tr2_inc_q4 /// Bus lane, income groups 1-4
	tr2_geo2 tr2_geo3 tr2_geo4 tr2_geo5 tr2_geo6 /// Bus lane, geography 2-6
	) order(toll_fam_mean /// Baseline toll effect
	tr1_inc_q1 tr1_inc_q2 tr1_inc_q3 tr1_inc_q4 /// Toll, income groups 1-4
	tr1_geo2 tr1_geo3 tr1_geo4 tr1_geo5 tr1_geo6 /// Toll, geography 2-6
	ptl_fam_km_mean /// Baseline bus lane effect
	tr2_inc_q1 tr2_inc_q2 tr2_inc_q3 tr2_inc_q4 /// Bus lane, income groups 1-4
	tr2_geo2 tr2_geo3 tr2_geo4 tr2_geo5 tr2_geo6 /// Bus lane, geography 2-6
	) se star(* 0.10 ** 0.05 *** 0.01)
	
* === The code below shows the right format of table D.2
* === but without stars to indicate significance, and without parantheses for SE

eststo clear
local outcome_vars BEV_fam_yes car_fam_count car_fam_yes ICE_fam_count
* For each outcome variable
foreach outcome in `outcome_vars' {
	* Load estimates
	estimates use "${ster}hetero_reg_O_`outcome'"
		* Store toll coefficients and standard errors as locals
		local toll1_`outcome' = _b[toll_fam_mean]
		local toll1_se_`outcome' = _se[toll_fam_mean]
		local toll2_`outcome' = _b[tr1_inc_q1]
		local toll2_se_`outcome' = _se[tr1_inc_q1]
		local toll3_`outcome' = _b[tr1_inc_q2]
		local toll3_se_`outcome' = _se[tr1_inc_q2]
		local toll4_`outcome' = _b[tr1_inc_q3]
		local toll4_se_`outcome' = _se[tr1_inc_q3]
		local toll5_`outcome' = _b[tr1_inc_q4]
		local toll5_se_`outcome' = _se[tr1_inc_q4]
		local toll6_`outcome' = _b[tr1_geo2]
		local toll6_se_`outcome' = _se[tr1_geo2]
		local toll7_`outcome' = _b[tr1_geo3]
		local toll7_se_`outcome' = _se[tr1_geo3]
		local toll8_`outcome' = _b[tr1_geo4]
		local toll8_se_`outcome' = _se[tr1_geo4]
		local toll9_`outcome' = _b[tr1_geo5]
		local toll9_se_`outcome' = _se[tr1_geo5]
		local toll10_`outcome' = _b[tr1_geo6]
		local toll10_se_`outcome' = _se[tr1_geo6]
		* Store bus lane coefficients and standard errors as locals
		local ptl1_`outcome' = _b[ptl_fam_km_mean]
		local ptl1_se_`outcome' = _se[ptl_fam_km_mean]
		local ptl2_`outcome' = _b[tr2_inc_q1]
		local ptl2_se_`outcome' = _se[tr2_inc_q1]
		local ptl3_`outcome' = _b[tr2_inc_q2]
		local ptl3_se_`outcome' = _se[tr2_inc_q2]
		local ptl4_`outcome' = _b[tr2_inc_q3]
		local ptl4_se_`outcome' = _se[tr2_inc_q3]
		local ptl5_`outcome' = _b[tr2_inc_q4]
		local ptl5_se_`outcome' = _se[tr2_inc_q4]
		local ptl6_`outcome' = _b[tr2_geo2]
		local ptl6_se_`outcome' = _se[tr2_geo2]
		local ptl7_`outcome' = _b[tr2_geo3]
		local ptl7_se_`outcome' = _se[tr2_geo3]
		local ptl8_`outcome' = _b[tr2_geo4]
		local ptl8_se_`outcome' = _se[tr2_geo4]
		local ptl9_`outcome' = _b[tr2_geo5]
		local ptl9_se_`outcome' = _se[tr2_geo5]
		local ptl10_`outcome' = _b[tr2_geo6]
		local ptl10_se_`outcome' = _se[tr2_geo6]		
	}

eststo clear
* For each row of the table
forvalues i = 1(1)10 {
	* Make a place to store the right coefficients
	eststo tr_`i' 
	* For each panel of the table (i.e. outcome)
	foreach outcome in `outcome_vars' {
		* Add scalars for toll coefficient and standard error
		estadd scalar `outcome'_toll_coef = `toll`i'_`outcome'',  replace
		estadd scalar `outcome'_toll_se   = `toll`i'_se_`outcome'',  replace
		* Add scalars for bus lane coefficient and standard error
		estadd scalar `outcome'_ptl_coef = `ptl`i'_`outcome'',  replace
		estadd scalar `outcome'_ptl_se   = `ptl`i'_se_`outcome'',  replace
		* Save
		eststo tr_`i'
	}
}	

* Print table (dimension 16*10)
esttab tr_*, keep("") ///
	stats( ///
		BEV_fam_yes_toll_coef BEV_fam_yes_toll_se BEV_fam_yes_ptl_coef BEV_fam_yes_ptl_se ///
		car_fam_count_toll_coef car_fam_count_toll_se car_fam_count_ptl_coef car_fam_count_ptl_se ///
		car_fam_yes_toll_coef car_fam_yes_toll_se car_fam_yes_ptl_coef car_fam_yes_ptl_se ///
		ICE_fam_count_toll_coef ICE_fam_count_toll_se ICE_fam_count_ptl_coef ICE_fam_count_ptl_se)
		
		